Tree Data

Code

devtools::install_github("dmurdoch/leaflet@crosstalk4")
library(tidyverse)
library(sf)
library(leaflet)
library(stringr)
library(knitr)
library(reactablefmtr)
library(plotly)
library(htmltools)
library(crosstalk)
library(patchwork)
library(ggtext)
library(DT)
library(showtext)

agg_df <- readRDS('../r_objects/agg_df.Rdata')
sa2_sf <- readRDS('../r_objects/sa2_sf.Rdata') %>%
  #rename(sa2_code_2021 = "SA2_CODE21") %>%
  st_set_geometry('geometry')

sal_sf <- readRDS('../r_objects/sal_sf.Rdata') %>%
  st_set_geometry('geometry')

source('../r_alt/functions.r')
source('../r/mapping_functions.r')
source('../themes/theme.r')

Tree coverage

Insert Lit Review here

Code

grouping = 'SAL_NAME21'
grouping_sf = sal_sf
pretty_group_name = 'Suburb'

trees_grouped <- coverage(group = grouping) %>%
  left_join(grouping_sf, by = (grouping)) %>%
  st_set_geometry('geometry') %>%
  dplyr::select(grouping, tree_percentage,distance, geometry) %>%
  rowwise() %>%
  mutate(tree_percentage = round(tree_percentage, 3)) %>%
  mutate(distance = distance / 1000) %>%
  as('Spatial')

#bar chart
sd_map = SharedData$new(trees_grouped)
sd_df = SharedData$new(as.data.frame(trees_grouped@data), group = sd_map$groupName())

  #mapCoverage(trees_per_sa2_cw, trees_per_sa2_cw$data()$tree_percentage , 'Greens', 'Tree Coverage (%)')

  treePal <- colorNumeric(palette = 'Greens', domain = as.data.frame(sd_df$data())$tree_percentage)
    
  leaflet(sd_map) %>%
    setView(lng = 144.963115, lat = -37.814175, zoom = 10) %>%
    addProviderTiles('CartoDB.Positron') %>%
    addPolygons(fillColor = ~treePal(tree_percentage),
                fillOpacity = 0.5,
                opacity = 0)  %>%
    addLegend(position = "bottomright",
              pal = treePal,
              values = sd_df$data()$tree_percentage,
              title = 'title')

Code

  filter_slider("distance", "Distance to CBD", sd_df, ~distance, min = 0, max = max(sd_df$data()$distance) + 1, width = '75%')

Distance to CBD

Code

  extra_cols = list()
  extra_cols[[grouping]] <- colDef(name = pretty_group_name, filterable = T)
         
  reactable(sd_df,
  defaultSorted = list('tree_percentage' = "desc"),
  pagination = TRUE,
  theme = sandstone(),
  columns = c(list(
    #grouping = colDef(filterable = TRUE),
    tree_percentage = colDef(
      name = 'Tree Coverage (%)',
      cell = data_bars(trees_grouped,
                       text_position = "above",
                       fill_color = c("#EDF8E9", "#006D2C"),
                       #fill_gradient = TRUE,
                       #fill_color_ref = viridi(5),
                       background = "lightgrey",
                       max_value = 100,
                       brighten_text = FALSE,
                       bar_height = 20
                       )
    ),
    distance = colDef(show = FALSE)
  ),extra_cols)
)

Coverage of streets, parks and residential areas

This section displays how much of each land use type is covered by trees. For example, a value of 60 for ‘Streets’ means that 60% of total street area is under canopy coverage. Accordingly, if ‘Parks’ is a low number, it means parks are not shaded, not that there are few parks in a region, it is just likely they are open spaces.

Code

street_trees_grouped <- coverage(exp = expression(zone_short == 'roads'), type = 'street', group = grouping)

reserve_trees_grouped <- coverage(exp = expression(feature_type == 'reserve'), type = 'reserve', group = grouping) 

resi_trees_grouped <- coverage(exp = expression(zoning_permits_housing == 'Housing permitted'), type = 'residential', group = grouping) 

#combine into one and clean
combined_df <- street_trees_grouped %>%
  left_join(reserve_trees_grouped, by = grouping) %>%
  left_join(resi_trees_grouped, by = grouping) %>%
  dplyr::select(c(grouping,'street_percentage', 'reserve_percentage', 'residential_percentage')) %>%
  left_join(grouping_sf %>% st_drop_geometry() %>% select(grouping, distance), by = grouping) %>%
  rowwise() %>%
  # mutate(distance = distance_map[as.character(as.name(grouping))]) %>%
  mutate(distance = distance / 1000) %>%
  mutate(across(where(is.numeric), ~ round(.x, 1))) %>%
  mutate(across(ends_with('percentage'), ~ ifelse(is.na(.x), 0, .x) )) %>%
  ungroup() %>%
  as.data.frame()

sd_combined = SharedData$new((combined_df))

Note: the table below is sorted by street tree coverage by default. Click other column names to sort them.

Note: There might be innaccuracies in the street coverage data when using Suburbs as grouping unit, as not all Statistical Area 1s fit uniquely within a suburb. As such, SA1s are assigned into suburbs based on the one they have the highest overlap with.

Code

filter_slider("distance", "Distance to CBD", sd_combined, ~distance, min = 0, max = max(sd_combined$data()$distance), width = '75%')

Distance to CBD

Code

#weird bubble plot  
reactable(sd_combined,
          defaultSorted = list('street_percentage' = "desc"),
          defaultColDef = colDef(
  cell = bubble_grid(combined_df, shape = 'circles', colors = c("#EDF8E9", "#006D2C")),
  align = 'center',
  vAlign = 'center'),
  columns = c(list(distance = colDef(show = F),
                 street_percentage = colDef(name = 'Streets'),
                 reserve_percentage = colDef(name = 'Parks'),
                 residential_percentage = colDef(name = 'Residential')),extra_cols),
  theme = sandstone(),
  defaultSortOrder = 'desc')

Interactions with Missing Middle

YM recently released Housing Targets, a data-driven model which allocates housing targets by council.

This page calculates how much tree cover might be lost in the upzoning of residential lots, and then calculating how to offset this loss by increasing coverage of street trees

Code

mm_lgas <- c('Brimbank', 'Merri-bek', 'Banyule', 'Darebin', 'Yarra', 'Moonee Valley', 'Manningham', 'Maribyrnong', 'Melbourne', 'Hobsons Bay', 'Port Phillip', 'Boroondara', 'Stonnington', 'Glen Eira', 'Bayside', 'Monash', 'Whitehorse', 'Maroondah', 'Manningham', 'Kingston')


years <- 10

targets <- read.csv('../council_targets.csv') [,1:2] %>%
  as.data.frame() %>%
  mutate(Target = sub(',', '', Target))

target_df <- data.frame()

for(lga in targets$LGA) {
  home_target <- (targets %>% filter(LGA == lga))$Target
 
  target_df <<- rbind(target_df, lga_targets(lga, as.numeric(home_target)*years))
}

colnames(target_df) = c('LGA', 'initial_coverage', 'final_coverage')

target_df <- target_df %>%
  rowwise() %>%
  mutate(percentage_point_change = as.numeric(final_coverage) - as.numeric(initial_coverage),
                                  percentage_change = ((as.numeric(final_coverage)/as.numeric(initial_coverage))-1)*100 ) %>%
  arrange(-percentage_change) %>%
  mutate(percentage_change = round(percentage_change, 3))

target_df %>% reactable(theme = sandstone())

Zoning

This section breaks down canopy coverage by zone type. As expected, low density residential and civic land (parks) top the list. It is tempting to note that Greenfield is below industrial, however, Greenfield by definition has not had enough time for planted trees to grow large enough for canopy cover.

Code

coverage(exp = expression(!(zone_short %in% c('roads','other'))),group = 'zone_short') %>%
  
  reactable(
    defaultSorted = list('tree_percentage' = "desc"),
    pagination = F,
    theme = sandstone(),
    columns = list(
    tree_percentage = colDef(
      name = 'Coverage %',
      cell = data_bars(.,
      text_position = "above",
                       fill_color = c("#EDF8E9", "#006D2C"),
                       #fill_gradient = TRUE,
                       #fill_color_ref = viridi(5),
                       background = "lightgrey",
                       max_value = 100,
                       brighten_text = FALSE,
                       bar_height = 20
    )),
    coverage = colDef(show = F),
    total_area = colDef(show = F),
    zone_short = colDef(name = 'Zone type')
  ))

Regressions

aaa todo: mix with Census data

Note: in the scatter plots below, suburbs with less than 500 lots suitable for residential development are excluded, alongside suburbs over 30km from the city.

Code

lot_sizes <- agg_df %>%
  filter(zoning_permits_housing == 'Housing permitted') %>%
  group_by(SAL_NAME21) %>%
  summarise(med_lot = median(lot_size, na.rm = TRUE), n = n()) %>%
  rowwise()

house_price_data = readxl::read_xls('../data/Median-House-Price-3rd-Quarter.xls', skip = 1)[-c(1:5),1:6] 

New names:
• `Jul - Sep 2022` -> `Jul - Sep 2022...2`
• `Apr - Jun 2023` -> `Apr - Jun 2023...5`
• `Jul - Sep 2023` -> `Jul - Sep 2023...6`
• `Jul - Sep 2023` -> `Jul - Sep 2023...7`
• `Jul - Sep 2022` -> `Jul - Sep 2022...9`
• `Apr - Jun 2023` -> `Apr - Jun 2023...10`

Code

house_price_data = house_price_data %>%
  rowwise() %>%
  mutate(avg = mean(as.numeric(
    c(`Jul - Sep 2022...2`,
      `Oct - Dec 2022`,
      `Jan - Mar 2023`,
      `Apr - Jun 2023...5`,
      `Jul - Sep 2023...6`)), na.rm = TRUE)) %>% 
  select(SUBURB, avg,) %>%
  rename(SAL_NAME21 = SUBURB) %>%
  mutate(SAL_NAME21 = str_to_title(SAL_NAME21)) %>%
  left_join(combined_df, by = grouping) %>%
  filter(!is.na(distance), !is.na(avg)) %>%
  rename(average_house_price = avg) %>%
  mutate(average_house_price = average_house_price / 10^6) %>%
  left_join(lot_sizes, by = 'SAL_NAME21') %>%
  adf() %>%
  filter(n > 500, distance < 30)

missing_suburbs = setdiff(combined_df$SAL_NAME21, house_price_data$SAL_NAME21)
#wtf

house_interactable = SharedData$new(house_price_data)

filter_slider("distance", "Distance to CBD", house_interactable, max = 30, ~distance, width = '75%')

Distance to CBD

Code

price_street_plot <- ggplot(house_interactable, mapping = aes(x = average_house_price, y = street_percentage, text = SAL_NAME21), axe) + xlab('Average House Price ($m)') + ylab('Street Tree Coverage (%)') + geom_point(color = "#006D2C") + theme_minimal()

price_residential_plot <- ggplot(house_interactable, mapping = aes(x = average_house_price, y = residential_percentage, text = SAL_NAME21), axe) + xlab('Average House Price ($m)') + ylab('Residential Area Tree Coverage (%)') + geom_point(color = "#006D2C") + theme_minimal()


#bscols(widths = c(4,4), ggplotly(price_street_plot), ggplotly(price_residential_plot))
subplot(ggplotly(price_street_plot),ggplotly(price_residential_plot), titleY = T, titleX = T)

Is this relationship driven by richer areas having larger residential lot sizes?

Code

lot_df = house_price_data %>% filter(med_lot < 2000)

lotsize_street_plot <- ggplot(lot_df, mapping = aes(x = med_lot, y = street_percentage, text = SAL_NAME21), axe) + xlab('Average Lot Size') + ylab('Street Tree Coverage (%)') + geom_point(color = "#006D2C") + theme_minimal()

lotsize_residential_plot <- ggplot(lot_df, mapping = aes(x = med_lot, y = residential_percentage, text = SAL_NAME21), axe) + xlab('Average Lot Size') + ylab('Residential Tree Coverage (%)') + geom_point(color = "#006D2C") + theme_minimal()


ggplotly(lotsize_street_plot)

Code

ggplotly(lotsize_residential_plot)

Code

#summary(lm(residential_percentage ~ med_lot, house_price_data)) 

#summary(lm(street_percentage ~ med_lot, house_price_data))